%dXt = sigma(t)^2dt 
%sigma(t)=OU
clear all
% index J for X
betaK=15; IIpart=sparse(index2(betaK));
% time T
T = 1/2;
% time mesh 
num=50; tspan = linspace(0,T,num+1);
%sigma
s=sigmahandle; a=1;b=1;m=1;sig0=1;
%options = odeset('AbsTol',1e-12,'RelTol',1e-6);
options = [];
row=1+2*betaK+betaK*(betaK-1)/2; X0=zeros(row,1); 
% solver of X
disp('solving X1')
tic
[~,X1] = ode45(@randucousqdt,tspan,X0,options,T,a,b,m,sig0,betaK,s,IIpart);

%dXt = sigma(t)dW 
%sigma(t)=OU

% index J for X
alphaK=20; %betaK=15; 
%sigma
sig=sigmahandle; %a=1;b=1;m=1;sig0=1;
s  =mihandle;
row=(betaK+1)*(alphaK+1); X0=zeros(row,1); 
% solver of X
disp('solving X2')
[time,X2] = ode45(@randucoudw,tspan,X0,options,T,a,b,m,sig0,alphaK,betaK,sig,s);


hold on
%%%the WCE price in matlab
price=zeros(51,1);
Npath=10000;
eta=randn(Npath,betaK);
xi=randn(Npath, alphaK);

for iter=1:51
asset=X2(iter,1)*ones(Npath,1);
for jj=1:betaK
    asset=asset+X2(iter,jj+1)*eta(:,jj);
end
for ii=1:alphaK
    asset=asset+X2(iter,ii*(betaK+1)+1)*xi(:,ii);
    for jj=1:betaK
        asset=asset+X2(iter,ii*(betaK+1)+1+jj)*xi(:,ii).*eta(:,jj);
    end
end
price(iter)=sum(abs(asset))/Npath;
end
toc
plot(time, price);



%%%the WCE price
price=zeros(51,1);
Npath=10000;
eta=randn(Npath,betaK);
xi=randn(Npath, alphaK);

for iter=1:51
asset=pricex2(eta,xi,X2(iter,:));
price(iter)=sum(abs(asset))/Npath;
end

plot(time, price,'r');

%%


hold on
%%%the WCE price in matlab
price=zeros(51,1);
Npath=100000;
eta=randn(Npath,betaK);
etasquare=1/sqrt(2)*(eta.^2-1);
etaeta=zeros(Npath, betaK*(betaK-1)/2);
counter=1;
for ii=1:betaK-1
    for jj=ii+1:betaK
        etaeta(:,counter)=eta(:,ii).*eta(:,jj);
        counter=counter+1;
    end
end

for iter=1:51
asset=zeros(Npath,1);
for ii=2:betaK+1
    asset=asset+X1(iter,ii)*eta(:,ii-1);
end
for ii=betaK+2:2*betaK+1
    asset=asset+X1(iter,ii)*etasquare(:,ii-betaK-1);
end
for ii=2*betaK+2:2*betaK+1+betaK*(betaK-1)/2
    asset=asset+X1(iter,ii)*etaeta(:,ii-2*betaK-1);
end
asset=asset+X1(iter,1);

price(iter)=sum(abs(asset))/Npath;
end
toc
plot(time, price);
%%%the WCE price in C
price=zeros(51,1);
Npath=10000;
eta=randn(Npath,betaK);
etasquare=1/sqrt(2)*(eta.^2-1);
etaeta=etaietaj(eta);
xi=randn(Npath, alphaK);

for iter=1:51
    asset=pricex1(eta,etasquare,etaeta,X1(iter,:));
    price(iter)=sum(abs(asset))/Npath;
end

toc
plot(time, price,'r');

%%
tic
Npath=5000;
mcX=zeros(5001,Npath);
dt=0.5/5000;
for pn=1:Npath
    disp(pn)
    samplepath=getapath(5001,'ou');
    for ii=2:5001
        mcX(ii,pn)=mcX(ii-1,pn)+samplepath(ii-1)^2*dt+samplepath(ii-1)*randn(1)*sqrt(dt);
    end
end
time=(0:50)/100;
tmp=sum(abs(mcX),2)/Npath;
mcprice=tmp(1:100:5001);
toc
plot(time,mcprice)

%%
%dXt = sigma(t)^2dt 
%sigma(t)=OU
clear all
% index J for X
betaK=15; IIpart=sparse(index2(betaK));
% time T
T = 1/2;
% time mesh 
num=50; tspan = linspace(0,T,num+1);
%sigma
s=sigmahandle; a=1;b=1;m=1;sig0=1;
%options = odeset('AbsTol',1e-12,'RelTol',1e-6);
options = [];
row=1+2*betaK+betaK*(betaK-1)/2; X0=zeros(row,1); 
% solver of X
disp('solving X1')
tic
[~,X1] = ode45(@randucousqdt,tspan,X0,options,T,a,b,m,sig0,betaK,s,IIpart);

%dXt = sigma(t)dW 
%sigma(t)=OU

% index J for X
alphaK=20; %betaK=15; 
%sigma
sig=sigmahandle; %a=1;b=1;m=1;sig0=1;
s  =mihandle;
row=(betaK+1)*(alphaK+1); X0=zeros(row,1); 
% solver of X
disp('solving X2')
[time,X2] = ode45(@randucoudw,tspan,X0,options,T,a,b,m,sig0,alphaK,betaK,sig,s);


%%%the WCE price
price=zeros(51,1);
Npath=10000;
eta=randn(Npath,betaK);
etasquare=1/sqrt(2)*(eta.^2-1);
etaeta=etaietaj(eta);
xi=randn(Npath, alphaK);

for iter=1:51
asset=pricex1(eta,etasquare,etaeta,X1(iter,:));
asset=asset+pricex2(eta,xi,X2(iter,:));
price(iter)=sum(abs(asset))/Npath;
end
toc
hold on
plot(time,price,'r')